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1 . Introduction 

In  the  past,  transonic  flows  have  been  calculated  with  the 

Euler  equations  and  with  time-dependent  finite-difference  methods.  For 

example,  Magnus  and  Yoshihara1  used  the  Lax-Wendroff  finite-difference 

schemes  with  additional  artificial  viscosity.  In  their  work,  the  shocks 

were  smeared,  and  small  grid  sizes  were  needed  to  capture  them.  Grossman 
2 

and  Moretti,  on  the  other  hand,  fitted  the  shock  following  a method 

3 

developed  by  Kentzer,  which  uses  the  compatibility  relations  along  the 
characteristics . 

For  many  cases  of  interest,  a potential  flow  model  is  adequate. 

4 

Murman  and  Cole  introduced  a type-dependent  finite-difference  scheme 
and  solved  the  transonic  small-disturbance  equation  by  relaxation  methods. 
Jameson,"*  using  "rotated  difference  schemes,"  extended  their  work  to  the 
full  potential  equation.  In  their  calculations,  shock  jump  conditions 
are  not  satisfied.  For  example,  mass  is  not  conserved  across  the  shock, 
and  the  strength  and  position  of  the  shock  are  usually  not  calculated 
correctly.  Later,  Murman^  introduced  a "fully  conservative  scheme"  to 
handle  this  problem  for  small-disturbance  calculations,  and  Jameson^ 
Introduced  a fully  conservative  scheme  for  the  full  potential  equation. 

In  these  later  calculations,  mass  is  conserved  globally,  supersonic  to 
subsonic  shocks  are  located  within  a few  grid  points,  and  supersonic  to 
supersonic  shocks  are  usually  smeared  over  more  grid  points.  Sharper 
shocks  can  be  obtained  only  by  grid  refinement. 

Q 

Hafez  and  Cheng  considered  shock  fitting  for  transonic  small- 

disturbance  calculations  by  using  type-dependent  finite-difference 

relaxation  methods.  In  their  work,  shock- jump  conditions  are  explicitly 

imposed,  and  mass  is  locally  conserved  across  a surface  of  discontinuity. 

Much  coarser  grids  may  then  be  used  for  the  calculations.  The  algorithm 

applies  for  supersonic  to  subsonic  shocks  as  well  as  for  supersonic  to 

supersonic  shocks.  Using  characteristic  compatibility  relations,  Yu  and 
9 

Seebass  studied  the  same  problem.  For  embedded  shocks,  they  used  a 
method  similar  to  that  of  Hafez  and  Cheng. 

Extension  of  the  method  of  Hafez  and  Cheng  to  the  full  potential 
calculation  is  straightforward  in  principle;  however,  in  practice,  the 
more  complicated  coordinate  systems  used  for  the  full  potential  calculations 
lead  to  cumbersome  interpolation  formulae. 
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In  this  paper  we  consider  an  alternative  procedure  for  fitting 
shock  waves  in  full  potential  calculations.  (The  application  of  the 
method  to  small  disturbances  is  discussed  in  the  Appendix.)  The  method 
is  based  on  an  equation  for  the  unsteady  shock  jump,  which  is  derived 
from  the  time-dependent  equation  describing  the  relaxation  algorithm. 

In  the  steady-state  limit,  shock  speed  vanishes,  and  the  steady  shock 
polar  is  retained.  In  this  way,  the  jump  conditions  are  imposed  iteratively 
and  in  a manner  consistent  with  the  relaxation  procedure  that  is  used 
everywhere  else  in  the  flow  field.  The  method  should  be  applicable  to 
rotated  difference  schemes  and  to  conservative  and  nonconservative 
differences . 

The  only  previous  work  reported  in  literature  on  shock-fitting 
methods  for  full  potential  calculations  is  the  work  of  Jones  and  South, ^ 
in  which  detached  bow  shock  waves  were  considered.  Jones  and  South  first 
used  a mapping  from  physical  space  to  a rectangular  computational  grid; 
then,  they  used  Newton's  method  to  determine  the  shock  shape.  The 
method  considered  in  this  paper  appears  to  be  simpler  and  is  not  restricted 
to  bow  shocks. 

In  the  following  we  first  discuss  unsteady  transonic  flow 
equations  and  their  weak  solutions.  Then,  the  time-dependent  equations 
describing  the  iterative  methods  are  then  developed.  A shock-fitting 
algorithm  is  described,  a solution  procedure  is  proposed,  and  some 
preliminary  numerical  results  are  given.  Finally,  in  the  Appendix  we 
compare  the  application  of  the  shock-fitting  algorithm  to  small-disturbance 
calculations  with  the  method  of  Hafez  and  Cheng  for  one-  and  two-dimensional 
numerical  examples. 
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2 . Unsteady  Full  Potential  Equations 

For  simplicity,  we  consider  here  the  Cartesian  coordinates  x , 
y and  t . The  velocity  potential  $ is  governed  by 


<J>  + 2u0  + 2v$  _ = (a2  - u2)c(i  - 2uv<J>  + (a2  - v2)<f> 

tt  xt  yt  xx  xy  yy 


where 


K + "2  + ’2 -1]  ■ 

n 


a) 


u = + cos  a , 


v = * + sin  a . 

y 


and  a is  the  angle  of  attack  of  the  Oncoming  flow  with  Mach  number 

H . 

00 

Equation  (1)  is  not  in  conservative  form.  Multiplying  by 

2 

p/a  , we  obtain  the  conservative  form 


-Pt  = (Pu)x  + (pv)y  , (2) 

where  p - ^1  - :L-~  — M2  (2cf>t  + u2  + v2  - 1)1  1/(Y  " 1)  . 

For  smooth  flows,  equations  (1)  and  (2)  are  equivalent.  Equation  (2), 
however,  admits  a weak  solution  with  mass  conservation  across  a discontinuity 
surface. 

2.1  Jump  Conditions 

The  jump  condition  admitted  by  the  weak  solution  is  given  by 

- St  i p D = G pu  D Sx  + Q pv  3 sy  , (3) 

where  | p J denotes  the  jump  is  p across  the  surface 

S(x,y, t)  = 0 . 

In  addition,  from  irrationality  conditions  we  have 
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(4) 


or  st  : sx  : sy " B *t  D : II  K D : 0 *y  D . 

which  implies 

I»tl  +g  I*J  +%  0 *y  D -o 


I <t>  D = o • 

(5) 

For  steady  states,  equations  (3) 

and  (5)  reduce  to 

avD 

Q pu  D - C pv  0 

■ 0 or 

0 p q„  0 ; 

(6) 

l U D f - + D V fl  = o or 

D qs  11  - o . 

(7) 

where  qn  is  the  relative  normal  velocity  to  the  shock  X - X° (y, t)  *»  0 , 
and  qg  is  the  tangential  velocity. 

For  transonic  flows,  equation  (6)  is  a good  approximation  of 
Rankine  Hugoniot  relations*,  the  difference  being  due  to  the  irrationality 
assumption. 


^tx  ™ *^xt  * 

and 

*ty  " *yt  * 

Thus, 

st  11 

♦x  0 " D D 

1 

♦y  11  - I *t  1 

The  exact  Frandtl  relation  (taking  into  consideration  entropy  variations)^ 
may  be  used  Instead  of  eq.  (6)  in  the  manner  suggested  by  Jones  and  South'; 
namely, 

Wd  K)u  * a*2  - % • 

a 

where  d and  u refer  to  downstream  and  upstream  conditions,  and  a is 
the  sonic  speed  of  sound. 
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2.2  Characteristics  and  Compatibility  Relations 

We  note  that  equation  (1)  or  (2)  is  always  hyperbolic.  ( t 
is  the  time-like  direction  for  both  subsonic  and  supersonic  flows.) 
Equation  (1)  can  be  written  in  the  following  form: 

+ 2q(J)  = (a2  - q2)di  + a2$  » 

Ttt  ^Yst  ’ rss  Ynn  * 


where 


and  q 
N = n , 
of  the 


4>  = — m (u^4>  + 2uv<)>  + v2<j>  ) , 

Tss  <2  ' Txx  Txy  Tyy 

<P  = (v24>  - 2uv<J>  + u2(J>  ) , 

nn  q2  Txx  xy  yy 


6 = — d)  + — rf)  . 

vst  q Txt  q vyt 


is  the  total  velocity.  Using  the  transformation 

■T  = t , equation  (8)  reduces  to:  (freezing  q , 

d)  . term) 
st 


S = s - qt  , 
the  coefficient 


|>TT  “ a (^gg  + 


kNN' 


(9) 


or,  with  the  transformation  X=x-ut,  Y = y - vt  , T = t , equation 
(1)  reduces  to 


<f> 


TT 


a ^XX  + ^YY^ 


(10) 


Equations  (1)  and  (4)  can  be  written  as  a system  of  first-order  equations; 
namely 


-uv 

0 

0 
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where 


W “ <J>t  , a2  = -^2  - (2w  + u2  + V2  - 1)  . (11) 

M 

00 

For  time-dependent  problems,  the  characteristics  and  the  compatibility 
relations  are  readily  obtained  from  this  system.  For  steady  problems, 
they  are  given  by 


where 


dv\  = _ A /djM  = _ B±  %/b2-AC 

du/  . C \dx/  , C 

ch  ch 


A 

B 

C 


(12) 
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3.  Shock-Fitting  Method 

3.1  Iterative  Solution  for  Steady  Equation 

It  is  well  known  that  iterative  line  relaxation  methods  can  be 
described  by  a time-dependent  equation.  For  example,  Jameson^  showed 
that  the  convergence  of  the  iterative  solution  of  the  full  potential 
equation  can  be  analyzed  by 


°*xt  + »yt  + ^yyt  + *t  = 


where  R^Ofr)  = (a2  - u2)^  - Zuvif)^  + (a2  - v2)<f> 


yy 


(13) 


Multiplying  equation  (9)  by  p/aZ  , we  have 


^t^yt  + V+^t  = RFcW)  ’ 


(14) 


where 


^(0)  “ (Pu)x  + (Pv) 


Since  we  are  interested  only  in  the  steady-state  solution,  the  coefficients 
ct  , 3 , Y and  6 may  be  selected  to  accelerate  the  convergence. 

The  jump  conditions  admitted  by  the  weak  solution  of  equation  (10) 
(freezing  the  coefficients  in  the  left-hand  side)  are  given  by 


st  C + % + Y<fyy  0 = C pu  D -sx  + A pv  D Sv  . 


(15) 


Let  the  shock  surface  S(x,y,t)  be  described  by  the  explicit  relation 
S = X - XD(y;t)  = 0 . 


Hence, 


fr  A p D = 0 pu  D - fr  B pv  I . 


(16) 


where 


p = -(ad)  + 3<p  + yd:  ) 

x t-Ty  ,Yyy' 


Equation  (16)  may  be  written  in  the  form 


AiE.+  B 3^- 

A 3t  + B 3y  “ C » ' 


(17) 
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where  A = |[  p j)  , 

b = n pv  d » and 

c “ C pu  U . 

i 

3.2  Shock-Fittine 

Shock-fitting  procedures  similar  to  methods  used  for  the 
unsteady  Euler  equations  may  be  applied  to  the  iterative  calculations 
once  the  artificial  time-dependent  equation  (17)  describing  the  develop- 
ment of  the  solution  through  iteration  is  recognized.  Since  we  are  only 
interested  in  the  steady-state  solution,  however,  many  simplifications 
may  be  made.  Two  cases  are  identified  here. 

3.2.1  Supersonic-Subsonic  Shocks 

Given  an  initial  condition  (an  initial  guess  of  the  shock 
location  XD(y;0)  and  a boundary  condition  (for  example,  3X°/3y  (0,t)  = 0 ) , 
equation  (17)  determines  the  new  shock  location  X^(y;t)  across  which 
<p  is  continuous.  Hence,  a Dirichlet  boundary  condition  is  imposed  on 
the  subsonic  (elliptic)  flow  downstream  of  the  shock  as  shown  in  Sketch  1. 


Y 


Marching  in  Time 


Sketch  1.  Supersonic-Subsonic  Shocks. 
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u 


3.2.2  Supersonic-Supersonic  Shocka 

In  general,  the  jump  conditions  and  the  characteristic  compati- 
bility relations  admitted  by  the  time-dependent  equation  uniquely  determine 
the  flow  just  downstream  of  the  shock  (<f>  and  <J>n)  . Similar  to  the 
supersonic-subsonic  shocks,  equation  (17)  is  used,  and  <|>  on  the  shock 
is  obtained  from  upstream  extrapolation.  However,  $ may  be  obtained 
by  solving  the  compatability  relations  along  the  downstream  characteristics 
together  with  the  shock  jump  conditions  (Equation  7) , simultaneously  to 
provide  u and  v downstream  of  the  shock.  These  are  needed  for  the 
supersonic  flow  calculations  since,  in  rotated  difference  schemes, 

‘Krv  * and  terms  contributing  to  t(>  are  backward  differenced. 

Hence,  the  values  of  <J>  at  two  upstream  grid  points  or  equivalently  the 

derivatives  <6  and  <p  are  needed.  We  encountered  some  difficulties 
x y 

in  implementing  this  scheme  for  supersonic-supersonic  shocks  in  general. 

The  details  of  a shock-fitting  algorithm  are  discussed  in  the  next  section. 
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4.  Numerical  Implementation  of  Shock-Fitting  Method 

A shock-fitting  algorithm  has  the  following  properties: 

(1)  An  initial  estimate  of  the  potential  is  obtained,  for 
example,  from  a nonconservative  solution. 

(2)  The  shock  waves  are  detected  and  used  as  an  initial 
estimate  of  shock  location. 

(3)  The  flow  field  is  computed  with  a relaxation  method,  in 
which  the  shock  is  a surface  of  discontinuity  fixed  in 
space  and  the  shock  jump  conditions  are  imposed  as  the 
boundary  conditions. 

(4)  The  shock  locations  are  updated,  based  on  the  latest 
information,  by  using  equation  (17). 

(5)  Steps  3 and  4 are  respected  until  convergence  is  achieved. 

Step  1 requires  no  elaboration  here.  For  step  2,  different  criteria  may 

be  used  to  detect  a shock  from  smooth  calculations.  Murman^  used  the 

maximum  slope  point  of  the  velocity  profile  in  the  shock  region.  South^ 

2 

used  the  minimum  Laplacian  (V  <j>)  for  full  potential  calculations.  For 
supersonic-subsonic  shocks,  the  shock  points  are  suitable  Indicators  of 
shock  locations.  Interpolation  may  be  used  between  mesh  points  to 
locate  the  sonic  velocity  for  a more  accurate  approximation.  Jones  and 
South^®  used  an  analytical  form  for  the  shock,  with  the  right  behavior 
in  the  far  field  (asymptotic  to  the  Mach  line)  and  near  the  axis  (normal 
shock).  To  determine  the  standoff  distance  of  a detached  bow  shock, 

Jones  and  South  used  the  point  with  the  same  Mach  number  as  the  one 
downstream  of  a normal  shock.  Here  we  use  a similar  procedure,  which  is 
discussed  in  Section  5.  For  Step  4,  the  new  shock  location  is  obtained 
as  a solution  of  equation  (17).  A Lax-Wendroff  finite-difference  scheme 
may  be  used  to  solve  this  equation.  To  avoid  a stability  restriction,  we 
also  tested  an  implicit  scheme  of  the  Crank-Nicholson  type.  A tridiagonal 
system  is  solved  each  time  the  shock  is  relocated.  The  finite-difference 
formula  are 
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Lax-Wendroff  Scheme: 


Here  we  neglected  the  variation  of  A,  B and  C with  time. 


Crank-Nicholson  Scheme: 


0.5 


Once  the  shock  Is  relocated,  the  continuity  of  $ is  imposed  as  a 

boundary  condition.  For  a detached  bow  shock,  (j)  is  set  equal  to  zero, 

the  free-stream  value.  For  embedded  shocks,  the  value  of  <f>  at  the 

shock  is  obtained  from  upstream  (supersonic)  conditions  by  extrapolation. 

In  the  supersonic-subsonic  case,  to  avoid  use  of  a special  unequal  mesh 

formula  at  the  grid  point  downstream  of  the  shock,  a Taylor  series 

expansion  is  used  to  obtain  a fictitious  value  of  $ at  the  first  mesh 

point  upstream  of  the  shock.  This  procedure  is  clearly  demonstrated  for 

the  one-dimensional  small-disturbance  case  in  the  appendix.  It  Is  used 

only  for  supersonic-subsonic  shocks  (i.e.,  if  the  flow  downstream  of  the 

shock  is  subsonic)  since  it  contradicts  the  rule  of  forbidden  signals  in 

the  case  of  downstream  supersonic  flow.  For  the  supersonic-supersonic 

shocks,  two  conditions  <f>  and  are  needed  to  continue  the  calculation 

n 

downstream  of  the  shock.  The  value  of  <)>  on  the  shock  is  obtained  from 
upstream  extrapolation;  can  be  obtained  from  studying  equations  7 

and  12  simultaneously.  We  have  not  succeeded  in  testing  the  scheme  for 
supersonic-supersonic  shocks. 
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5.  Calculated  Examples 

We  have  obtained  preliminary  results  for  the  flow  past  a blunt 

12 

body  by  using  the  method  of  South  and  Jameson  with  the  above-described 

14 

shock-fitting  algorithm.*  We  use  the  Keller-South  (RAXBOD)  program  to 

obtain  an  Initial  estimate  for  $ everywhere  in  the  flow  field.  The 

standoff  distance  Is  calculated  as  discussed  earlier,  and  hyperbola  is 

used  as  an  initial  estimate  of  the  shock  shape. 

The  shock  is  located  by  using  equation  (11).  Properties  ahead 

of  the  shock  are  calculated  by  extrapolation  from  upstream  conditions; 

In  this  case,  c|)  = 0 ahead  of  the  shock.  For  the  supersonic-supersonic 

part  of  the  shock,  we  also  need  the  properties  behind  the  shock.  Equations 

(7)  and  (12)  provide  $ and  <p  (hence  $ ) just  downstream  of  the 

x y q 

shock.  In  this  example,  we  assume  that  ^ (i.e.,  the  derivative  of  <J> 

along  the  normal  computational  coordinate  direction)  just  downstream  of 
the  shock  is  negligible.  More  accurate  treatment  for  the  supersonlc- 
aupersonlc  shocks  are  meeded.  Figure  1 shows  the  results  for  the  pressure 
coefficient  on  the  body  calculated  by  the  RAXBOD  program  for  a sphere  of 
radius  one  in  a flow  with  = 1.2  . The  physical  plane  is  mapped  to  a 
rectangular  grid  (21  x 21  grid  points),  where 

AY 

n — - , a = 1.3  , A o 1.25  . 

<1  - Y)“ 

Here,  n is  the  physical  coordinate  normal  to  the  body,  and  Y is  the 
computational  coordinate,  which  varies  from  zero  at  the  body  to  one  at 
infinity.  The  tangential  coordinate  is  stretched  by  a quadratic  transfor- 
mation between  the  physical  arc  length  and  the  computational  coordinate 

X . 


The  same  problem  was  solved  by  Jones  and  South,  who  used  mapping 
techniques,  and  by  Hsieh1^  (hemisphere  cylinder),  who  used  the  time- 
dependent  Euler  equation. 
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Figure  1.  Cp  on  a Sphere  at  M = 1.2  (RAX BOD  Solution). 


Figure  2 shows  the  initial  guess  of  the  shock  location.  The 
movement  of  the  shock  (standoff  distance)  with  iteration  is  plotted  in 
Figure  3,  and  the  final  result  is  shown  in  Figure  4. 

Shock  Line 


-1.8864  3.1974  8.2813  13.3651  18.4490  23.5328 

Figure  2.  Initial  Shock  Location  Detected  From  RAXBOD  Solution. 
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Figure  4.  Shock  Location  After  Shock  Fitting. 
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It  should  be  mentioned  here  that  our  results  Cor  the  case  of  a 
sphere  at  = 1.2  (based  on  the  full  potential  calculations)  indicate 
that  the  bow  shock  standoff  distance  differs  by  20  percent  from  a time- 
dependent  solution  of  Euler  Equations  given  in  Reference  13.  These 
discrepancies  cannot  be  traced  to  the  difference  between  the  Euler 
Equation  and  the  full  potential  equation  since  in  this  range  they  should 
agree. 

It  could  be  argued  that  the  discrepancies  are  due  to  the 
supersonic/  supersonic  part  of  the  shock  where  the  shock  was  not  really 
fitted  in  our  calculation.  This  is  unlikely,  however,  since  the  initial 
guess  is  selected  using  Moeckel's  approximation.  We  have  not  resolved 
this  point  and  further  investigation  is  needed  to  resolve  the  differences. 

Calculation  of  embedded  shocks  have  been  tested  also.  Flows 
around  a sphere  were  computed  with  and  without  shock  fitting.  In  Figure 
5,  results  of  the  RAXBOD  program  are  plotted.  For  the  same  mesh,  results 
from  the  SAXBOD  program  are  plotted  in  Figure  6.  Notice  that  the  overall 
pressure  distribution  does  not  seem  to  look  different,  however,  the  drag 
coefficient  using  Trapezoidal  and  Simpson's  Rule  are  changed  by  about 
20  percent.  The  same  case  has  been  calculated  after  a mesh  refinement. 

The  results  of  the  RAXBOD  program  are  shown  in  Figure  7,  while  the  results 
of  the  SAXBOD  program  are  shown  in  Figures  8 and  9. 


Cp 
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Figure  5.  Cp  on  a Sphere  at  M = 1.2  (RAXBOD  Solution). 
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Figure  6.  Cp  on  a Sphere  at  M = 1.2  (SAXBOD  Solution). 
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Figure  7.  Cp  on  a Sphere  at  M = 1 .2. 


20 


N> 


1.1704 

1.1688 

1.1620 

1.1468 

1.1199 

1.0778 

1.0172 

.9355 

.8306 

.7018 

.9492 

.3746 

.1810 

-.0274 

-.2454 

-.4682 

-.0908 

-.9069 

-1.1106 

-1.2969 

-1.4638 

-1.6097 

-1.7332 

-1.8357 

-1.9187 

-1.9846 

-2.0513 

-1.7092 

-.1060 

.7495 

.7657 

.8315 

.9099 

.9840 

1.0463 

1.0945 

1.1288 

1.1510 

1.1635 

1.1691 

1.1704 


Figure  8.  Cp  on  a Sphere  at  M = 1.2. 
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Figure  9.  Cp  on  a Sphere  at 


M = 1.2. 
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In  the  first  case,  the  iterations  were  stopped  when  the  maximum 

-4 

difference  between  successive  iterates  is  less  than  10  while  more 
iterations  were  allowed  (to  10  for  the  results  shown  in  Figure  5.  The 
maximum  residual  (in  the  difference  equation  or  the  shock  polar)  is 

_3 

less  than  10  . The  difference  in  the  pressure  distribution  is  indis- 

tinguishable, but  the  drag  coefficient  shows  2 percent  changes.  Comparing 
RAXBOD  and  SAXBOD  results,  the  drag  coefficients  are  increased  by  mesh 
refinement  when  shock  fitting  is  used,  and  decreased  in  the  case  of 
RAXBOD  (nonconservative  calculations) . It  would  be  interesting  to  run 
the  same  case  using  a conservative  code.  The  computer  listing  of  SAXBOD 
for  embedded  shocks  are  attached. 


\ 


23 


AEDC-TR-78-37 


6.  Conclusions 

A shock-fitting  algorithm  for  full  potential  type-dependent 
finite-difference  relaxation  calculations  Is  described,  and  preliminary 
results  are  presented.  The  method  Is  based  on  a time-dependent  equation 
describing  the  line  relaxation  method.  Written  in  conservative  form, 
its  weak  solution  admits  an  unsteady  jump  condition  in  terms  of  the 
shock  speed.  We  use  this  equation  to  update  the  shock  during  the  iterative 
calculations.  Numerical  implementations  of  the  method  are'  discussed, 
and. the  flow  around  a sphere  is  calculated.  We  have  hot  succeeded  in 
satisfactorily  implementing  the  shock  fitting  algorithm  for  the  supersonic- 
supersonic  shock  and  more  accurate  treatment  is  needed.  We  also  test 
this  method  for  small-disturbance  calculations. 
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Appendix 

i 

In  this  section,  we  consider  the  application  of  the  proposed  algorithm 
to  the  small-disturbance  equation.  The  unsteady  small-disturbance 
equations  may  be  written  as  (0  = a = 1): 

+ 2atiXt = 


d>  = d> 
xt  *tx 


^yt  = ^ty  ’ 


Hence,  the  jump  conditions  admitted  by  the  weak  solution  are 


(*£  HU  +2aM»  0*j)st" 

0(1  - *£>♦,  - **Cy  + Dl&jjl  sx  + C4>yD  sy 
st  O^B  - 0*tD  sx 
st  hj  - DM  sx  . 

For  0 = 0,  a = 1 , S(x,y,t)  = X - XD(y,t)  = 0 . We  have5" 

*£  Ir  - - < » - £ - <y  + »«^>+ 

' lr  0 1 - 0 *t ) . 

‘ # I ♦,  I - - I ♦,  i if  • 


Equation  (A. 3)  may  be  written  in  the  equivalent  (quasisteady)  form: 

0 (i  - *1)  - a&J  - [*  I . 0 


(A.l) 


(A. 2) 


(A. 3) 
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For  the  steady  state,  equation  (A. 3)  reduces  to  the  shock  polar 


- <«  - «l)  - (Y  + 1)  - (§) 


2 

shock 


Equation  (A.l)  may  he  written  as  a system  of  first-order  equations 
(w  = (j>t  , u - (j>x  , v = 4>y) : 

ewt  = [(1  - m£)  u - H(y  + 1)  M2u2]  x 

+ v - a2M2w 
y <»x 


(A. 4) 


vt  = wy  • (A. 5) 

Notice  that,  unlike  the  full  potential  case,  the  system  of  first-order 
equations  is  hyperbolic  only  if 

“V  + 1 - M2  > 0 where  1 - M2  = 1 - M2  - (y  + l)M2d> 

OO  ' w OO 


The  characteristics  (and  the  compatibility  relations)  of  the 
unsteady  equations  are  obtained  from  this  system  (A. 5).  The  steady- 
state  relations  are: 

(0)  = (£)  = (Y  + 1)  M2  u - (1  - M2)  . (A.6) 

' char.  ' ' char. 

We  note  that  equation  (A. 4)  simplifies  to* 


Equation  (A.6)  is  consistent  with  the  weak  shock  approximation:  The  shock 
bisects  the  characteristics;  in  other  words. 
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(A. 7) 


The  development  of  solution  during  Iteration  may  be  described  by 


a<j,xt  + ^yyt  + = 

f(i  - m2)<)>  - + (4>  ) 

L “ x 2 »yxj  x Yy  i 


Hence,  the  shock  Is  updated  from  the  equation  (for  y = 0) 

a D . D 
A 3t“  + B 37"  = C * 


where 


A “ ct  , B = 


3x 

3y  * 


C = - <(1  - M2)  - (Y  + 1)M2®  > . 


One-Dimensional  Example 

Consider  the  following  one-dimensional  model  equation  for  the 
unsteady  problem 


(A. 8) 


(A. 9) 


(4>x>t  + " 0 * 0 1 * < * (A. 10) 

♦ (0,t)  = 0 , 6 (0,t)  = <f> 

<KA,t)  = 0 


The  jump  condition  is 


3X°  , 

■57—  B <{>  + <}> 

3t  *L  % 

We  see  that  the  speed  of  the  shock  is  the  average  of  the  upstream  and 
downstream  velocities. 


(A. 11) 


The  governing  equation  is  ut  + 2uux  = 0 ; hence 

(£)  " 2u  and  (£)  „ = <2uu  + 2ud>  ■ 
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For  Z - 2 and  <f>  = 1 , the  steady-state  solution  is 

*L 

4>  = x o<x<i  , 

(A. 12) 

<f>  = 2 - X 1 < X < 2 , 

In  Sketch  A.l  we  see  that  if  we  locate  the  shock  incorrectly  between 
zero  and  one  the  shock  will  move  towards  one.  If  the  shock  is  placed 
at  one,  the  shock  speed  is  zero,  and  the  shock  settles  there. 


(I)  (II) 


Sketch  A.1 


To  check  these  ideas,  we  compared  this  method  to  that  of  Hafez  and 

g 

Cheng  . Hafez  and  Cheng  considered  the  shock  an  internal  boundary  with 
a derivative  (Neuman)  boundary  condition.  Assuming  that  the  initial 
guess  of  the  shock  location  is  at  I,  we  can  solve  the  flow  downstream  of 
the  shock  with  a derivative  boundary  condition  derived  from  the  steady 
jump  condition;  in  other  words, 


(A. 13) 


The  new  position  of  the  shock  is  determined  by  finding  the  intersection 
of  the  two  $ solutions  downstream  and  upstream  of  I,  as  demonstrated 
in  in  Sketch  A. 2. 
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The  continuity  of  <j>  across  the  shock  (tangential  velocities  are 

the  same)  and  the  mass  conservation  condition  switch  roles  in  the  two 

8 

methods.  Hafez  and  Cheng  used  the  mass  conservation  condition  to 
calculate  the  flow  downstream  of  the  shock  and  the  continuity  of  <f>  to 
determine  the  new  shock  location.  In  the  present  method  we  use  the 
continuity  of  <j>  to  calculate  the  flow  downstream  of  the  shock,  and  the 
shock  moves  to  satisfy  the  mass  conservation  condition.  Hafez  and  Cheng 
show  that  the  results  obtained  with  their  method  are  Identical  to  those 
obtained  with  Hurman's  fully  conservative  shock-point  operator  for  the 
one-dimensional  case.  On  the  other  hand,  the  present  method  uses  the 
Dlrichlet  boundary  condition  to  solve  the  flow  downstream  of  the  shock. 
Unlike  the  nonconservative  calculations,  mass  is  conserved  by  moving  the 
shock  to  the  correct  place. 

The  finite-difference  implementation  of  the  Hafez  and  Cheng  procedure 
for  the  one-dimensional  problem  considered  here  is 


AX 


0 


(A.  14) 


where 


*1+1  " *i 
AX 

*i-l  ~ *i-2 
AX 


» 


while  the  corresponding  formula  of  the  present  method  is  (Sketch  A. 4) 


where 


2 2 
u+  “ u- 


AX 


- 0 , 


(A. 15) 


u+  = 


u 


*i+i  : *i 

AX 

*i  - *d 
AX 


♦d  ■ *1-1  + ^ - *1-1' 
t *1-1  ' *i-2  *1+1  ’ *i.\ 

V AX  " AX  f ' 
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and 


dXD  ^1+1  ~ *i  ^i-1  ~ ^i-2 
dt  AX  AX 


which  can  be  written  in  the  form: 


u+  - u- 

AX 


(A. 16) 


Notice  the  left-hand  side  is  the  nonconservative  difference  scheme  (elliptic 
operator)  and  the  right-hand  side  is  the  correction  (source)  term  needed 
to  conserve  mass. 


♦ 


X 
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Two-Dimensional  Examples 


Supersonic-Subsonic  Shocks. 

Q 

In  the  method  of  Hafez  and  Cheng  , an  algebraic  relation  derived 
from  a finite-difference  approximation  of  the  shock  polar  replaces  the 
finite-difference  approximation  of  the  differential  equation  downstream 
of  the  shock;  namely. 


<1  - - <Y  + l)Mjf>x> 


bj  Q^yD 

Ax  + Ay 


0 


and  <J>  are  used 
x y 

m tne  u jj  and  < > quantities.  The  shock  is  located  at  the  inter- 
section of  the  $ surfaces  extrapolated  from  upstream  and  downstream 
calculations  Invoking  the  continuity  of  <j>  . 

In  the  present  method,  an  elliptic  operator  with  a Correction 
(source)  term  is  used  at  the  grid  points  just  downstream  of  the  shock. 
Fictitious  points  upstream  of  the  shock  similar  to  those  in  the  one- 
dimensional case  (as  shown  in  Sketch  A. 5)  are  needed.  The  values  of  <p 
at  these  points  ((f)^  , are  obtained  by  extrapolation  from  the 

previous  downstream  solution  and  from  the  values  of  <j>  at  the  shock. 

The  latter  are  obtained  by  extrapolation  of  the  upstream  solution  to  the 
shock  position.  The  shock  is  then  relocated  by  means  of  equation  (A. 9). 
Note  that  here  we  need  the  extrapolation  procedures  only  to  calculate  a 
correction  term.  (An  approximate  location  of  the  shock,  for  example, 
the  middle  of  the  mesh,  may  be  used  without  great  loss  of  accuracy. ) 
Also,  a locally  normal  shock-fitting  approximation  may  be  assumed  (in 
other, words,  we  may  use  the  same  formula  as  in  the  one-dimensional  case, 

plus  a centered  difference  approximation  for  the  <Ji  term).  As  a 
4 yy 

matter  of  fact,  Murman's  fully  conservative  scheme  may  be  written  as  an 

elliptic  operator  with  a correction  term,  as  follows: 


where 


fl  = tL  M. 
0 " Ax  3y 


D 


kL 

Ax 


<0 


and  where  appropriate  one-sided  derivatives  for 


rr  n 


Murman’s  shock-point  operator  is  a special  case  of  this  relation: 

(3  = 1 and  0 $y  II  My  = WyyJc.D^ 
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K - + *yy]c.D.  ' ‘K  ' • 

All  of  these  variants  have  been  tested  by  modifying  Hurman's  small- 
disturbance  program.  As  shown  in  figure  10,  we  have  calculated  the  flows 
around  a 6%  parabolic  arc  airfoil  for  a subsonic  free  stream  (M  = 0.9) 
by  using 

(i)  Nonconservative  scheme 
(ii)  Hurman's  fully  conservative  scheme 
(ill)  Fully  conservative  scheme  in  the  form  of  an  elliptic  operator 
with  a correction  (source)  term 
(iv)  Shock  fitting  (present  method). 

For  these  calculations,  we  used  a 60  x 60  grid  (Ax  = Ay  = 0.05),  and  the 
far  field  was  set  equal  to  zero.  The  convergence  limit  1 - <Jin|  < 10  ^ 
was  usually  reached  in  50  iterations. 


Shock 


i 

i 

X 

Sketch  A.5 


Supersonic-Supersonic  Shocks. 

g 

For  supersonic-supersonic  shocks,  Hafez  and  Cheng  used  equation  (A. 7) 
to  locate  the  shock  in  a step-by-step  procedure.  With  the  present 
method,  we  solve  the  problem  Iteratively  where  the  shock  moves  according 
to  equation  (A.  9).  The  value  of  <j>  at  the  shock  is  used  as  a boundary 
condition  for  the  unknown  on  a vertical  line  between  the  shock  location 
and  the  body  (axis).  Since  background  differences  are  used,  values  along 
the  i-1  and  i-2  lines  are  needed.  Extrapolation  in  the  Y direction 
is  used  to  obtain  the  values  of  4 for  those  points  upstream  of  the 
shock.  As  shown  in  figure  11,  we  have  calculated  the  flows  around  a 6% 
parabolic  arc  airfoil  for  a supersonic  stream  (M^  = 1.15)  by  using 
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Figure  10.  Supersonic-Subsonic  Shocks. 


SONIC  LINE  CdCHiMilolE' 

_J Y i.Sl'L.K .. 


J l7 

3? 

n 

34 

35 

L- 

espouse 

._"0;7 1. 

J 6 

' )2 

ii 

”34" 

~3S" 

2 

.05606 

■ i U U 7 t 

.<>'j  7; " 

15 

32 

31 

34 

35 

3 

.10000 

. Jlttfl  l 

. ttn fa 

)4 

32 

33 

34 

35 

36 

4 

. . 1 5900  . 

. -a  1 i •>!* 

, 744411 

17 

3? 

33 

34 

35 

J4 

5 

.70060 

• 1127  1 

.75011 

1? 

32 

3 1 

34 

35 

36 

6 

.75600 

. l4-»0  7 

. 76696 

11 

32 

33 

34 

35 

36 

7 

. 1(1  006 

« >ri£»6r'_ 

.X52-J. 

'io 

3?” 

"Yf 

■ 34" 

P— 35  76’ 

H 

• 35000 

•404PU 

• t Jll'l  f> 

9 

32 

33 

34 

35 

36 

9 

.4000(1 

• *»  ? 0 6 '/ 

. 7 

6 

32 

33 

34 

35 

36 

37 

.13 

.4.560  0 

■«» 

*745-4  l 

7 

32 

3 i 

36 

17 

1) 

Isuooo 

.*•4700 

• /5tf40 

6 

32 

33 

3* 

35 

36 

37 

12 

.55000 

.*•  74f>0 

_ S 

.12 

33 

3<» 

35 

36 

37 

Jl_ 

9^95qFi 

"4 

32 

“33 

3<* 

*”35' 

36 

37 

14 

■65C00 

.13  056 

3 

32 

33 

34 

35 

36 

17 

1 5 

• 7(j000 

."39S2 

.71459 

?. 

12 

33 

34 

35 

36 

_37_  | 

ln 

. ,15994  .. 

__  tS44!'0. 

.45604  . 

1 

j2 

14“ 

35 

16 

37  ' 

1 t 

.<’0(100 

.00225 

.67423 

Figure  10(a). 

Nonconservative  Calculations. 
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Figure  10(c).  Modified  Shock  Point  Operator  Results 

(Elliptic  Operator  With  a Correction  Term). 
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Figure  10(d).  Shock-Fitting  Results. 
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Figure  11.  Supersonic-Supersonic  Shocks 
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Figure  11{a).  Nonconservative  Calculations  (Did  Not  Converge  in  200  Iterations). 
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Figure  11(d).  Shock-Fitting  Results. 
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Figure  11(c).  Initial  Shock  Location  Predicted  From  Nonconservative  Calculations. 
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(i)  Nonconservative  scheme 
(ii)  Fully  conservative  scheme 
(ill)  Shock  fitting  (present  method); 

For  the  calculations  shown  in  figure  11,  we  used  a 60  x 30  grid 
(Ax  = 0.05  , Ay  = 0.1). 
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